This notebook is the relative cumulative frequency analysis of metabolic cage data.

Data was generated previously when I ran PCA, see notebook .

The idea behind this was to apply a method to study large metabolic datasets quantitatively. Based on the Gordon paper, we can use relative cumulative frequency to find frequency of data binned by VO2 or RQ reading levels, then calculate EC50. We can use EC50 values for every mouse in a simple t-test or ANCOVA to compare groups. In other words, this can be thought of as a kind of weighted average.

  • Halatchev IG, O’Donnell D, Hibberd MC, Gordon JI. Applying indirect open-circuit calorimetry to study energy expenditure in gnotobiotic mice harboring different human gut microbial communities. Microbiome. 2019;7(1):158. Published 2019 Dec 12. doi:10.1186/s40168-019-0769-4
  • Riachi M, Himms-Hagen J, Harper ME. Percent relative cumulative frequency analysis in indirect calorimetry: application to studies of transgenic mice. Can J Physiol Pharmacol. 2004;82(12):1075‐1083. doi:10.1139/y04-117

There are some assumptions that need to be considered. First are the flaws in the dataset — it is difficult to draw any conclusions when the datasets are flawed to begin with. There are 12 SPF mice, 5 of which are WT and 7 of which are L-Bmal1-KO. There are 6 GF mice, of which mouse #4 was contaminated and cannot be used, and of the remaining 5, 2 are WT and 3 are L-Bmal1-KO. The mix of genotypes prevents any far-reaching conclusions. Furthermore, due to limitations in the tools I am using (there needs to be the same number of samples for both groups), I am stuck to using only 5 SPF mice, to compare with 5 GF mice. Overcoming this requirement (i.e. having unequal groups) will be a good idea in the future.

Thus, this notebook serves as a proof-of-concept, that relative cumulative frequency can be applied, EC50 calculated, and statistics run.

Relative cumulative frequency

See the paper for more information . This is a method where data is first ordered by increasing value, then checked for increases based on some increment, e.g. for every increment with an increase, a frequency is recorded. Finally, frequency is cumulated. This is all accomplished by the function rcf. In order to run this, however, the range of data needs to be checked. find_range can be run to check the range on a vector, and will return a recommended range and increment value to submit to rcf.

EC50

By fitting frequencies to intervals, a sigmoidal regression is found. This is dependent on a normal distribution — that is, the mouse should have less frequencies towards the extremes, and the most expression in the middle. For RQ and VO2, this is true, the mouse will not spend as much time on the extremes of the values. For other data types, the nature and distribution of the data should be considered before using in this analysis.

t-test and ANCOVA

For the Welch t-test, I run an F-test first (var.test), which should be not significant.

dataset

Overview

These datasets are from RMR data collected for mice after introduction to the metabolic cages (days 0-6, see image).

overview

SPF

Start of
Met Cage run		 Metabolic
Cage #		 Ear Tag DOB Gender Genotype Body Mass on
BMR day (g)
3/10/2020 1 1555 11/18/19 M WT 26.4
3/10/2020 2 1692 12/06/19 M WT 25.6
3/10/2020 3 1693 12/06/19 M WT 24.7
3/10/2020 4 1695 12/06/19 M WT 25.6
3/10/2020 5 1699 12/06/19 M WT 28.7
3/10/2020 6 1556 11/18/19 M L-Bmal1-KO 26.2
3/10/2020 9 1557 11/18/19 M L-Bmal1-KO 27.1
3/10/2020 10 1558 11/18/19 M L-Bmal1-KO 26.9
3/10/2020 11 1559 11/18/19 M L-Bmal1-KO 24.7
3/10/2020 12 1691 12/06/19 M L-Bmal1-KO 28
3/10/2020 13 1694 12/06/19 M L-Bmal1-KO 26.4
3/10/2020 14 1700 12/06/19 M L-Bmal1-KO 27.7

GF

Cage Ear Tag Genotype GF condition, day 6 GF condition, day 10
1 5650 L-Bmal1-KO GF contaminated
2 976 WT GF GF
3 977 WT GF contaminated
4 990 L-Bmal1-KO likely contaminated likely contaminated
5 979 L-Bmal1-KO GF contaminated
6 978 L-Bmal1-KO GF contaminated

For the first half of the experiment, only cage 4 had contamination. For the second leg of the experiment, cages 1,3,5 and 6 were contaminated, and it is likely cage 4 was also contaminated.

See the 16S rDNA gel below.

gel

libraries

reqpkg <- c("docstring","magrittr","dplyr","ggplot2","ggpubr","reshape2","tibble","scales","purrr","nplr", "car","compute.es","effects","multcomp","pastecs","WRS2")
for (i in reqpkg) {
    print(i)
    print(packageVersion(i))
    suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "docstring"
## [1] '1.0.0'
## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "reshape2"
## [1] '1.4.3'
## [1] "tibble"
## [1] '2.1.3'
## [1] "scales"
## [1] '1.1.0'
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'
theme_set(theme_pubr())

functions

signif.floor <- function(x, n){
  pow <- floor( log10( abs(x) ) ) + 1 - n
  y <- floor(x / 10 ^ pow) * 10^pow
  # handle the x = 0 case
  y[x==0] <- 0
  y
}

signif.ceiling <- function(x, n){
  pow <- floor( log10( abs(x) ) ) + 1 - n
  y <- ceiling(x / 10 ^ pow) * 10^pow
  # handle the x = 0 case
  y[x==0] <- 0
  y
}

find_range <- function(data, channel, segments=10) {
  data <- data %>% dplyr::select(contains(channel)) %>% melt()
  r <- range(data$value)
  cut <- max(data$value)-min(data$value)
  out <- list("range"=r, "cut"=cut/segments)
  return(out)
}

rcf <- function(data, channel, id, min, max, b=0.1) {
  #' Relative cumulative frequency function
  #'
  #' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
  #'
  #' @param data numeric vector
  #' @param channel character. Channel to label the rcf output.
  #' @param b double. Number to cut vector.
  #' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
  
  breaks <- seq(min, max, by=b)
  # breaks <- seq(signif.floor(min(data), 2), signif.ceiling(max(data), 2), by=b)
  duration.cut <- cut(data, breaks, right = F)
  duration.freq <- table(duration.cut)
  duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
  duration.cumrelfreq = duration.cumfreq/length(data)
  out <- duration.cumrelfreq %>% as.data.frame()
  out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
  out$id <- rep(id, nrow(out))
  return(out)
}

EC50 <- function(x) {
  out <- lapply(x, function(y) {
    nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
  }) %>% 
    sapply(function(z) {
      getEstimates(z, 0.5)$x
    })
  out %<>% set_rownames(NULL)
  return(out)
}

lm_eqn <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                   list(a = format(unname(coef(m)[1]), digits = 2),
                        b = format(unname(coef(m)[2]), digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));
}

lm_eqn_str <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
  return(eq)
}

remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

ANCOVA <- function(df) {
  # Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
  
  cat("linear fit of dependent variable with covariates\n")
  p1 <- ggplot(df, aes(df[[2]], df[[3]])) +
    geom_smooth(method = "lm", se = FALSE) +
    geom_point() +
    # geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    ggtitle("Linear model - DV and Covar") +
    xlab("DV") + ylab("Covar")
  print(p1)
  
  m <- lm(df[[3]] ~ df[[2]], df)
  
  summary(m) %>% print()
  lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
  
  plot1 <- ggplot(df, aes(x=df[[1]], y=df[[2]])) +
    geom_boxplot() +
    # ggtitle("DV vs IV") +
    xlab("IV") + ylab("DV")
  plot2 <- ggplot(df, aes(x=df[[1]], y=df[[3]])) +
    geom_boxplot() +
    # ggtitle("Covar vs IV") +
    xlab("IV") + ylab("Covar")
  p2 <- ggarrange(plot1, plot2,
                  labels = c("DV vs IV", "Covar vs IV"),
                  ncol = 2)
  # grid.arrange(plot1, plot2, ncol=2)
  print(p2)
  
  cat("\n_____________________\n")
  cat("Some stats on the variables\n")
  cat("DV - IV:\n")
  by(df[[2]], df[[1]], stat.desc) %>% print()
  cat("Covar - IV:\n")
  by(df[[3]], df[[1]], stat.desc) %>% print()
  cat("levene's test\n")
  leveneTest(df[[2]], df[[1]], center = median) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of DV with IV\n")
  model <- aov(df[[2]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of covariate with IV, independent of DV\n")  
  model <- aov(df[[3]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of DV with IV, normalized to covariate\n")
  c <- df[[2]]/df[[3]]
  model <- aov(c ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANCOVA\n")
  model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
  summary(model) %>% print()
  out <- Anova(model, type= "III")
  print(out)
  return(out)
}

check_ANCOVA <- function(res) {
  if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
  } else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
    print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
  } else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
  }
}

import data

SPF_dark <- read.csv2("data/R/SPF_dark.csv", sep = ",") %>% 
  set_colnames(paste0("SPF_",colnames(.)))
SPF_light <- read.csv2("data/R/SPF_light.csv", sep = ",") %>% 
  set_colnames(paste0("SPF_",colnames(.)))
SPF_covar <- read.csv2("data/R/SPF_covar.csv", sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

GF_dark <- read.csv2("data/R/GF_dark.csv", sep = ",") %>% 
  dplyr::select(-contains(c("_M_4","_R_4","_Mnz_4"))) %>% 
  set_colnames(paste0("GF_",colnames(.))) 
GF_light <- read.csv2("data/R/GF_light.csv", sep = ",") %>% 
  dplyr::select(-contains(c("_M_4","_R_4","_Mnz_4"))) %>% 
  set_colnames(paste0("GF_",colnames(.)))
GF_covar <- read.csv2("data/R/GF_covar.csv", sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

light <- SPF_light %>% 
  bind_cols(GF_light) %>% 
  dplyr::select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

dark <- SPF_dark %>% 
  bind_cols(GF_dark) %>% 
  dplyr::select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

total <- bind_rows(light, dark)

VO2

Hide

Select a tab above to view

Diurnal cycle

total.VO2 <- total %>% dplyr::select(contains("VO2"))
total.VO2.SPF <- total.VO2 %>% dplyr::select(contains("SPF"))
total.VO2.GF <- total.VO2 %>% dplyr::select(contains("GF"))

# Only taking last 5 KO mice
for (name in colnames(total.VO2.SPF)) {
  if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
    total.VO2.SPF %<>% dplyr::select(-name)
  }
}

rcf.SPF <- lapply(total.VO2.SPF, function(x) rcf(x, "VO2", "SPF", 1, 2.9, 0.1))
rcf.GF <- lapply(total.VO2.GF, function(x) rcf(x, "VO2", "GF", 1, 2.9, 0.1))

data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)

Models <- lapply(data, function(tmp){
  nplr(tmp$breaks, tmp$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice VO2 in complete diurnal cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 1.68268673163679"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 1.54667034693964"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.4691, num df = 4, denom df = 4, p-value = 0.7184
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.1529618 14.1102621
## sample estimates:
## ratio of variances 
##           1.469126
t.test(t[,1], t[,2])
## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = 2.1814, df = 7.7213, p-value = 0.06193
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.008805302  0.284815805
## sample estimates:
## mean of x mean of y 
##  1.680105  1.542100
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7915 -0.5634 -0.0086  0.3746  3.6986 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  26.7421     7.8325   3.414  0.00917 **
## df[[2]]       0.4766     4.8496   0.098  0.92413   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.733 on 8 degrees of freedom
## Multiple R-squared:  0.001206,   Adjusted R-squared:  -0.1236 
## F-statistic: 0.009659 on 1 and 8 DF,  p-value: 0.9241
## 
## y = 27 + 0.48x, r^2 = 0.00121

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  5.000000000  0.000000000  0.000000000  1.394711940  1.611752055  0.217040115 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  7.710499431  1.593015270  1.542099886  0.040261519  0.111783896  0.008104949 
##      std.dev     coef.var 
##  0.090027493  0.058379806 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000   1.57230931   1.86329413   0.29098482 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   8.40052569   1.65195795   1.68010514   0.04879998   0.13549046   0.01190719 
##      std.dev     coef.var 
##   0.10912007   0.06494836 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  27.10000000  31.20000000   4.10000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 141.40000000  27.50000000  28.28000000   0.75193085   2.08769472   2.82700000 
##      std.dev     coef.var 
##   1.68136849   0.05945433 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 133.70000000  26.90000000  26.74000000   0.58360946   1.62035962   1.70300000 
##      std.dev     coef.var 
##   1.30499042   0.04880293 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0129 0.9123
##        8               
## _____________________
## ANOVA of DV with IV
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.04761 0.04761   4.758 0.0607 .
## Residuals    8 0.08005 0.01001                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1  5.929   5.929   2.618  0.144
## Residuals    8 18.120   2.265               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value  Pr(>F)   
## df[[1]]      1 1.698e-04 1.698e-04   15.31 0.00446 **
## Residuals    8 8.868e-05 1.109e-05                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.04761 0.04761   5.491 0.0516 .
## df[[3]]      1 0.01935 0.01935   2.232 0.1788  
## Residuals    7 0.06070 0.00867                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.008611  1  0.9931 0.35217  
## df[[1]]     0.066812  1  7.7053 0.02746 *
## df[[3]]     0.019352  1  2.2319 0.17883  
## Residuals   0.060696  7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0274620158378014"

Light cycle

light.VO2 <- light %>% dplyr::select(contains("VO2"))
light.VO2.SPF <- light.VO2 %>% dplyr::select(contains("SPF"))
light.VO2.GF <- light.VO2 %>% dplyr::select(contains("GF"))

# Only taking last 5 KO mice
for (name in colnames(light.VO2.SPF)) {
  if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
    light.VO2.SPF %<>% dplyr::select(-name)
  }
}

rcf.SPF <- lapply(light.VO2.SPF, function(x) rcf(x, "VO2", "SPF", 1, 2.9, 0.1))
rcf.GF <- lapply(light.VO2.GF, function(x) rcf(x, "VO2", "GF", 1, 2.9, 0.1))

data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)

Models <- lapply(data, function(tmp){
  nplr(tmp$breaks, tmp$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice VO2 in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 1.4791478793743"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 1.35106555611578"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.6667, num df = 4, denom df = 4, p-value = 0.6328
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.1735299 16.0076157
## sample estimates:
## ratio of variances 
##           1.666674
t.test(t[,1], t[,2])
## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = 2.1057, df = 7.5294, p-value = 0.07048
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01294961  0.25464656
## sample estimates:
## mean of x mean of y 
##  1.472360  1.351512
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8906 -0.7058 -0.2072  0.6747  3.4101 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   31.480      7.540   4.175   0.0031 **
## df[[2]]       -2.812      5.327  -0.528   0.6119   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.704 on 8 degrees of freedom
## Multiple R-squared:  0.03366,    Adjusted R-squared:  -0.08713 
## F-statistic: 0.2787 on 1 and 8 DF,  p-value: 0.6119
## 
## y = 31 + -2.8x, r^2 = 0.0337

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  5.000000000  0.000000000  0.000000000  1.241538599  1.450636943  0.209098344 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  6.757558245  1.375201548  1.351511649  0.035144752  0.097577474  0.006175768 
##      std.dev     coef.var 
##  0.078586054  0.058146783 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000   1.36405635   1.61172432   0.24766797 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   7.36180063   1.48720545   1.47236013   0.04537177   0.12597224   0.01029299 
##      std.dev     coef.var 
##   0.10145437   0.06890595 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  27.10000000  31.20000000   4.10000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 141.40000000  27.50000000  28.28000000   0.75193085   2.08769472   2.82700000 
##      std.dev     coef.var 
##   1.68136849   0.05945433 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 133.70000000  26.90000000  26.74000000   0.58360946   1.62035962   1.70300000 
##      std.dev     coef.var 
##   1.30499042   0.04880293 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   0.344 0.5737
##        8               
## _____________________
## ANOVA of DV with IV
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.03651 0.03651   4.434 0.0683 .
## Residuals    8 0.06588 0.00823                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1  5.929   5.929   2.618  0.144
## Residuals    8 18.120   2.265               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0001274 1.274e-04   9.084 0.0167 *
## Residuals    8 0.0001122 1.402e-05                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.03651 0.03651   3.985 0.0861 .
## df[[3]]      1 0.00174 0.00174   0.189 0.6765  
## Residuals    7 0.06414 0.00916                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.026050  1  2.8431 0.13563  
## df[[1]]     0.034800  1  3.7980 0.09232 .
## df[[3]]     0.001736  1  0.1895 0.67647  
## Residuals   0.064139  7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)

Dark cycle

dark.VO2 <- dark %>% dplyr::select(contains("VO2"))
dark.VO2.SPF <- dark.VO2 %>% dplyr::select(contains("SPF"))
dark.VO2.GF <- dark.VO2 %>% dplyr::select(contains("GF"))

# Only taking last 5 KO mice
for (name in colnames(dark.VO2.SPF)) {
  if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
    dark.VO2.SPF %<>% dplyr::select(-name)
  }
}

rcf.SPF <- lapply(dark.VO2.SPF, function(x) rcf(x, "VO2", "SPF", 1, 2.9, 0.1))
rcf.GF <- lapply(dark.VO2.GF, function(x) rcf(x, "VO2", "GF", 1, 2.9, 0.1))

data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)

Models <- lapply(data, function(tmp){
  nplr(tmp$breaks, tmp$VO2, silent = TRUE, useLog = F)
  })

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice VO2 in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
abline(h=0.5)

print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 1.90041457533305"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 1.73951498204456"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.1638, num df = 4, denom df = 4, p-value = 0.8867
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.1211708 11.1776467
## sample estimates:
## ratio of variances 
##           1.163789
t.test(t[,1], t[,2])
## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = 2.1973, df = 7.9544, p-value = 0.05943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.009573524  0.388474314
## sample estimates:
## mean of x mean of y 
##  1.919810  1.730359
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7701 -0.5941 -0.0510  0.4014  3.6662 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  25.8554     6.4773   3.992    0.004 **
## df[[2]]       0.9066     3.5364   0.256    0.804   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.727 on 8 degrees of freedom
## Multiple R-squared:  0.008148,   Adjusted R-squared:  -0.1158 
## F-statistic: 0.06572 on 1 and 8 DF,  p-value: 0.8041
## 
## y = 26 + 0.91x, r^2 = 0.00815

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000   1.55076022   1.85131885   0.30055863 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   8.65179656   1.80215775   1.73035931   0.05861441   0.16273970   0.01717825 
##      std.dev     coef.var 
##   0.13106581   0.07574485 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000   1.78102178   2.14647377   0.36545199 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   9.59904853   1.90591446   1.91980971   0.06323268   0.17556205   0.01999186 
##      std.dev     coef.var 
##   0.14139256   0.07364926 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  27.10000000  31.20000000   4.10000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 141.40000000  27.50000000  28.28000000   0.75193085   2.08769472   2.82700000 
##      std.dev     coef.var 
##   1.68136849   0.05945433 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 133.70000000  26.90000000  26.74000000   0.58360946   1.62035962   1.70300000 
##      std.dev     coef.var 
##   1.30499042   0.04880293 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   2e-04 0.9902
##        8               
## _____________________
## ANOVA of DV with IV
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.08973 0.08973   4.828 0.0592 .
## Residuals    8 0.14868 0.01859                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1  5.929   5.929   2.618  0.144
## Residuals    8 18.120   2.265               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)   
## df[[1]]      1 0.0002812 2.812e-04   17.06 0.0033 **
## Residuals    8 0.0001318 1.648e-05                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.08973 0.08973   6.323 0.0401 *
## df[[3]]      1 0.04934 0.04934   3.477 0.1045  
## Residuals    7 0.09934 0.01419                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.001463  1  0.1031 0.75754  
## df[[1]]     0.137125  1  9.6624 0.01712 *
## df[[3]]     0.049339  1  3.4766 0.10451  
## Residuals   0.099341  7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0171186600523453"

RQ

Hide

Select a tab above to view

Diurnal cycle

total.RQ <- total %>% dplyr::select(contains("RQ"))
total.RQ.SPF <- total.RQ %>% dplyr::select(contains("SPF"))
total.RQ.GF <- total.RQ %>% dplyr::select(contains("GF"))

# Only taking last 5 KO mice
for (name in colnames(total.RQ.SPF)) {
  if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
    total.RQ.SPF %<>% dplyr::select(-name)
  }
}

rcf.SPF <- lapply(total.RQ.SPF, function(x) rcf(x, "RQ", "SPF", 0.6, 1.1, 0.05))
rcf.GF <- lapply(total.RQ.GF, function(x) rcf(x, "RQ", "GF", 0.6, 1.1, 0.05))

data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice RQ in complete diurnal cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 0.840277055484309"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 0.846032542612996"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 0.27142, num df = 4, denom df = 4, p-value = 0.2345
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.02826007 2.60690625
## sample estimates:
## ratio of variances 
##          0.2714247
t.test(t[,1], t[,2])
## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = 0.047799, df = 6.0224, p-value = 0.9634
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02793109  0.02904508
## sample estimates:
## mean of x mean of y 
## 0.8444767 0.8439197
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.53056 -0.25067 -0.04439  0.39201  1.76874 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    81.49      20.60   3.955   0.0042 **
## df[[2]]       -63.95      24.40  -2.621   0.0306 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.272 on 8 degrees of freedom
## Multiple R-squared:  0.4619, Adjusted R-squared:  0.3947 
## F-statistic: 6.868 on 1 and 8 DF,  p-value: 0.03062
## 
## y = 81 + -64x, r^2 = 0.462

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 5.0000000000 0.0000000000 0.0000000000 0.8098007597 0.8735776070 0.0637768474 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 4.2195986627 0.8471038023 0.8439197325 0.0103345386 0.0286932792 0.0005340134 
##      std.dev     coef.var 
## 0.0231087309 0.0273826171 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 5.0000000000 0.0000000000 0.0000000000 0.8274327561 0.8574982941 0.0300655380 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 4.2223836496 0.8485680324 0.8444767299 0.0053841326 0.0149487485 0.0001449444 
##      std.dev     coef.var 
## 0.0120392864 0.0142565046 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  27.10000000  31.20000000   4.10000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 141.40000000  27.50000000  28.28000000   0.75193085   2.08769472   2.82700000 
##      std.dev     coef.var 
##   1.68136849   0.05945433 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 133.70000000  26.90000000  26.74000000   0.58360946   1.62035962   1.70300000 
##      std.dev     coef.var 
##   1.30499042   0.04880293 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.6428 0.4459
##        8               
## _____________________
## ANOVA of DV with IV
##             Df    Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.0000008 0.0000008   0.002  0.963
## Residuals    8 0.0027158 0.0003395               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1  5.929   5.929   2.618  0.144
## Residuals    8 18.120   2.265               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 7.240e-06 7.236e-06   1.579  0.244
## Residuals    8 3.667e-05 4.584e-06               
## _____________________
## ANCOVA
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0000008 0.0000008   0.005 0.9457  
## df[[3]]      1 0.0016246 0.0016246  10.422 0.0145 *
## Residuals    7 0.0010912 0.0001559                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##                Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.0278748  1 178.8168 3.065e-06 ***
## df[[1]]     0.0003705  1   2.3769   0.16704    
## df[[3]]     0.0016246  1  10.4221   0.01449 *  
## Residuals   0.0010912  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with Covar, p = 0.0144853254058274"

Light cycle

light.RQ <- light %>% dplyr::select(contains("RQ"))
light.RQ.SPF <- light.RQ %>% dplyr::select(contains("SPF"))
light.RQ.GF <- light.RQ %>% dplyr::select(contains("GF"))

# Only taking last 5 KO mice
for (name in colnames(light.RQ.SPF)) {
  if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
    light.RQ.SPF %<>% dplyr::select(-name)
  }
}

rcf.SPF <- lapply(light.RQ.SPF, function(x) rcf(x, "RQ", "SPF", 0.6, 1.1, 0.05))
rcf.GF <- lapply(light.RQ.GF, function(x) rcf(x, "RQ", "GF", 0.6, 1.1, 0.05))

data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice RQ in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 0.786092284744787"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 0.789466997689693"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 0.19478, num df = 4, denom df = 4, p-value = 0.1421
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.020280 1.870769
## sample estimates:
## ratio of variances 
##          0.1947798
t.test(t[,1], t[,2])
## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -0.10582, df = 5.5013, p-value = 0.9195
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04945103  0.04543731
## sample estimates:
## mean of x mean of y 
## 0.7849262 0.7869331
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.89818 -0.60505 -0.06311  1.03107  1.81545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    52.56      13.40   3.923   0.0044 **
## df[[2]]       -31.88      17.04  -1.871   0.0983 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.446 on 8 degrees of freedom
## Multiple R-squared:  0.3044, Adjusted R-squared:  0.2174 
## F-statistic: 3.501 on 1 and 8 DF,  p-value: 0.09825
## 
## y = 53 + -32x, r^2 = 0.304

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  5.000000000  0.000000000  0.000000000  0.727127187  0.833409714  0.106282527 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  3.934665528  0.791403096  0.786933106  0.017350154  0.048171749  0.001505139 
##      std.dev     coef.var 
##  0.038796123  0.049300408 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 5.0000000000 0.0000000000 0.0000000000 0.7571566743 0.8014503624 0.0442936881 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 3.9246312361 0.7869148767 0.7849262472 0.0076572940 0.0212600564 0.0002931708 
##      std.dev     coef.var 
## 0.0171222298 0.0218138072 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  27.10000000  31.20000000   4.10000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 141.40000000  27.50000000  28.28000000   0.75193085   2.08769472   2.82700000 
##      std.dev     coef.var 
##   1.68136849   0.05945433 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 133.70000000  26.90000000  26.74000000   0.58360946   1.62035962   1.70300000 
##      std.dev     coef.var 
##   1.30499042   0.04880293 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.1621 0.3125
##        8               
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000010 0.0000101   0.011  0.918
## Residuals    8 0.007193 0.0008992               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1  5.929   5.929   2.618  0.144
## Residuals    8 18.120   2.265               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 5.390e-06 5.393e-06   1.038  0.338
## Residuals    8 4.154e-05 5.193e-06               
## _____________________
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.000010 0.0000101   0.017 0.8992  
## df[[3]]      1 0.003109 0.0031092   5.329 0.0543 .
## Residuals    7 0.004084 0.0005834                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##                Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.0302126  1 51.7844 0.0001781 ***
## df[[1]]     0.0009266  1  1.5883 0.2479580    
## df[[3]]     0.0031092  1  5.3292 0.0543059 .  
## Residuals   0.0040840  7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)

Dark cycle

dark.RQ <- dark %>% dplyr::select(contains("RQ"))
dark.RQ.SPF <- dark.RQ %>% dplyr::select(contains("SPF"))
dark.RQ.GF <- dark.RQ %>% dplyr::select(contains("GF"))

# Only taking last 5 KO mice
for (name in colnames(dark.RQ.SPF)) {
  if (as.numeric(strsplit(name, "_")[[1]][4]) < 10) {
    dark.RQ.SPF %<>% dplyr::select(-name)
  }
}

rcf.SPF <- lapply(dark.RQ.SPF, function(x) rcf(x, "RQ", "SPF", 0.6, 1.1, 0.05))
rcf.GF <- lapply(dark.RQ.GF, function(x) rcf(x, "RQ", "GF", 0.6, 1.1, 0.05))

data <- bind_rows(bind_rows(rcf.SPF), bind_rows(rcf.GF)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of SPF vs GF mice RQ in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("SPF EC50: ", getEstimates(Models$SPF, 0.5)$x))
## [1] "SPF EC50: 0.905230988953924"
print(paste0("GF EC50: ", getEstimates(Models$GF, 0.5)$x))
## [1] "GF EC50: 0.884114999600456"
t <- list(rcf.SPF, rcf.GF) %>% sapply(EC50) %>% set_rownames(NULL)
var.test(t[,1], t[,2])
## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.3281, num df = 4, denom df = 4, p-value = 0.79
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.1382837 12.7562514
## sample estimates:
## ratio of variances 
##           1.328149
t.test(t[,1], t[,2])
## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = 2.2069, df = 7.8442, p-value = 0.05902
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.001406623  0.059374039
## sample estimates:
## mean of x mean of y 
## 0.9079685 0.8789848
df <- data.frame(meta=c(rep("SPF", nrow(t)), rep("GF", nrow(t))), ch=c(t[,1],t[,2]), covar=c(SPF_covar[c(8,9,10,11,12),], GF_covar[-4,]))
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2052 -0.4948 -0.2523  0.2112  1.6742 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    77.19      11.13   6.933  0.00012 ***
## df[[2]]       -55.61      12.46  -4.464  0.00210 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.928 on 8 degrees of freedom
## Multiple R-squared:  0.7135, Adjusted R-squared:  0.6777 
## F-statistic: 19.93 on 1 and 8 DF,  p-value: 0.0021
## 
## y = 77 + -56x, r^2 = 0.714

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 5.0000000000 0.0000000000 0.0000000000 0.8572243491 0.8966261382 0.0394017891 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 4.3949240720 0.8868342514 0.8789848144 0.0086072838 0.0238976511 0.0003704267 
##      std.dev     coef.var 
## 0.0192464717 0.0218962506 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  5.000000000  0.000000000  0.000000000  0.878609328  0.940941881  0.062332554 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  4.539842610  0.907323949  0.907968522  0.009919496  0.027540936  0.000491982 
##      std.dev     coef.var 
##  0.022180667  0.024428894 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  27.10000000  31.20000000   4.10000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 141.40000000  27.50000000  28.28000000   0.75193085   2.08769472   2.82700000 
##      std.dev     coef.var 
##   1.68136849   0.05945433 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   5.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 133.70000000  26.90000000  26.74000000   0.58360946   1.62035962   1.70300000 
##      std.dev     coef.var 
##   1.30499042   0.04880293 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0244 0.8798
##        8               
## _____________________
## ANOVA of DV with IV
##             Df  Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.00210 0.0021001    4.87 0.0584 .
## Residuals    8 0.00345 0.0004312                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1  5.929   5.929   2.618  0.144
## Residuals    8 18.120   2.265               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 2.047e-05 2.047e-05   3.581 0.0951 .
## Residuals    8 4.574e-05 5.717e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.002100 0.0021001   11.24 0.0122 *
## df[[3]]      1 0.002142 0.0021421   11.47 0.0117 *
## Residuals    7 0.001308 0.0001868                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.031750  1 169.9716 3.638e-06 ***
## df[[1]]     0.000282  1   1.5107   0.25874    
## df[[3]]     0.002142  1  11.4673   0.01166 *  
## Residuals   0.001308  7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with Covar, p = 0.0116573479895595"

session info

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] WRS2_1.0-0       pastecs_1.3.21   multcomp_1.4-12  TH.data_1.0-10  
##  [5] MASS_7.3-51.5    survival_3.1-11  mvtnorm_1.1-0    effects_4.1-4   
##  [9] compute.es_0.2-4 car_3.0-7        carData_3.0-3    nplr_0.1-7      
## [13] purrr_0.3.3      scales_1.1.0     tibble_2.1.3     reshape2_1.4.3  
## [17] ggpubr_0.2.5     ggplot2_3.3.0    dplyr_0.8.5      magrittr_1.5    
## [21] docstring_1.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] splines_3.6.3     assertthat_0.2.1  cellranger_1.1.0  yaml_2.2.1       
##  [5] pillar_1.4.3      lattice_0.20-40   glue_1.3.2        digest_0.6.25    
##  [9] ggsignif_0.6.0    minqa_1.2.4       colorspace_1.4-1  sandwich_2.5-1   
## [13] cowplot_1.0.0     htmltools_0.4.0   Matrix_1.2-18     survey_3.37      
## [17] plyr_1.8.6        pkgconfig_2.0.3   haven_2.2.0       openxlsx_4.1.4   
## [21] rio_0.5.16        lme4_1.1-21       mgcv_1.8-31       farver_2.0.3     
## [25] withr_2.1.2       nnet_7.3-13       klippy_0.0.0.9500 cli_2.0.2        
## [29] crayon_1.3.4      readxl_1.3.1      evaluate_0.14     fansi_0.4.1      
## [33] nlme_3.1-144      forcats_0.5.0     xml2_1.2.5        foreign_0.8-75   
## [37] tools_3.6.3       data.table_1.12.8 hms_0.5.3         mitools_2.4      
## [41] lifecycle_0.2.0   stringr_1.4.0     munsell_0.5.0     zip_2.0.4        
## [45] compiler_3.6.3    rlang_0.4.5       grid_3.6.3        nloptr_1.2.2.1   
## [49] labeling_0.3      rmarkdown_2.1     boot_1.3-24       gtable_0.3.0     
## [53] codetools_0.2-16  abind_1.4-5       DBI_1.1.0         reshape_0.8.8    
## [57] roxygen2_7.1.0    curl_4.3          R6_2.4.1          zoo_1.8-7        
## [61] knitr_1.28        mc2d_0.1-18       stringi_1.4.6     Rcpp_1.0.4.6     
## [65] vctrs_0.2.4       tidyselect_1.0.0  xfun_0.12
 

By Sumeed Yoyo Manzoor

smanzoor@uchicago.edu



Back to top